A synthesis of dryland restoration techniques.

Purpose

To quantitatively examine the efficacy of vegetation restoration in drylands globally.

Questions

1. What is the global extent of research that directly examined restoration of drylands?
   
2. What were the common measures?
   
3. Is the restoration of vegetation a common and primary focus?
   
4. How frequently does the restoration measure outcomes beyond the focal species?
   
5. What were the primary restoration goals as reported by primary authors?
   
6. How much variation was there in the techniques tested and how long were experiments monitored and tested?
   
7. How relatively effective were the techniques?

Step 1. Search

Study-level viz to document patterns in exclusions primarily and the relatie frequenices, at the study level, of major categories of evidence.

#study data####
library(tidyverse)
studies <- read_csv("data/studies.csv")
studies

## # A tibble: 278 x 18
##       ID title technique data  region exclude rationale observations
##    <dbl> <chr> <chr>     <chr> <chr>  <chr>   <chr>     <chr>       
##  1   152 Shor~ seeding,~ expe~ Africa no      <NA>      <NA>        
##  2   180 Rest~ chemical~ App.~ Africa no      <NA>      <NA>        
##  3   229 Infl~ soil see~ fiel~ Africa no      <NA>      <NA>        
##  4   230 Acti~ planting  fiel~ Africa no      <NA>      <NA>        
##  5   255 The ~ grazing ~ fiel~ Africa no      <NA>      <NA>        
##  6   262 Reve~ seeding,~ eper~ Africa no      <NA>      <NA>        
##  7   263 The ~ phytogen~ fiel~ Africa no      <NA>      <NA>        
##  8   264 Eval~ seeding,~ fiel~ Africa no      <NA>      <NA>        
##  9   271 Patc~ natural ~ fiel~ Africa no      <NA>      <NA>        
## 10     4 Fact~ natural ~ App.~ Africa no      <NA>      <NA>        
## # ... with 268 more rows, and 10 more variables: disturbance <chr>,
## #   system <chr>, goal <chr>, intervention <chr>, paradigm <chr>,
## #   grazing <chr>, hypothesis <chr>, soil <chr>, benchmark <chr>,
## #   notes <chr>

#quick look at rationale needed
exclusions <- studies %>%
  filter(exclude == "yes")

#quick look at studies with paradigms
evidence <- studies %>%
  filter(exclude == "no")

#library(skimr)
#skim(evidence)

#study-level viz#####
#exclusions
#ggplot(exclusions, aes(rationale, fill = region)) +
 # geom_bar() +
 # coord_flip() +
 # labs(x = "rational for exclusion", y = "frequency") +
 # scale_fill_brewer(palette = "Paired")

ggplot(evidence, aes(disturbance, fill = paradigm)) +
  geom_bar(na.rm = TRUE) +
  coord_flip() +
  scale_fill_brewer(palette = "Paired") +
  labs(y = "frequency")

#ggplot(evidence, aes(region, fill = paradigm)) +
 # geom_bar(na.rm = TRUE) +
 # coord_flip() +
 # scale_fill_brewer(palette = "Paired") +
 # labs(y = "frequency")

#ggplot(evidence, aes(data, fill = paradigm)) +
 # geom_bar(na.rm = TRUE) +
  #coord_flip() +
  #scale_fill_brewer(palette = "Paired") +
  #labs(y = "frequency")

#ggplot(evidence, aes(system, fill = paradigm)) +
 # geom_bar(na.rm = TRUE) +
 # coord_flip() +
 # scale_fill_brewer(palette = "Paired") +
 # labs(y = "frequency")

#ggplot(evidence, aes(goal, fill = paradigm)) +
  #geom_bar(na.rm = TRUE) +
  #coord_flip() +
  #scale_fill_brewer(palette = "Paired") +
  #labs(x = "outcome", y = "frequency")

#step 1 models####
#paradigm
derived.evidence <- evidence %>%
  group_by(technique, data, region, disturbance, goal, paradigm) %>% summarise(n = n())

#active-passive split
m <- glm(n~paradigm, family = poisson, derived.evidence)
anova(m, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                       167     9.9147         
## paradigm  1 0.045115       166     9.8696   0.8318

#region
m1 <- glm(n~paradigm*region, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##                 Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                              167     9.9147         
## paradigm         1 0.045115       166     9.8696   0.8318
## region           6 0.301367       160     9.5682   0.9995
## paradigm:region  6 0.213627       154     9.3546   0.9998

#outcome
m2 <- glm(n~paradigm*goal, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m2, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                            167     9.9147         
## paradigm       1 0.045115       166     9.8696   0.8318
## goal           6 0.240941       160     9.6287   0.9997
## paradigm:goal  4 0.301480       156     9.3272   0.9897

#even split between active and passive evidence by all key categories

Step 2. Sort

A summary of sort process using PRISMA.

library(PRISMAstatement)
prisma(found = 1504,
       found_other = 5,
       no_dupes = 1039, 
       screened = 1039, 
       screen_exclusions = 861, 
       full_text = 178,
       full_text_exclusions = 101, 
       qualitative = 77, 
       quantitative = 66,
       width = 800, height = 800)

Step 3. Synthesize

Check data and calculate necessary measures.

#all data
#data <- read_csv("data/data.csv")
#data <- data %>%
  #mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/n.t*mean.t^2) + (sd.c^2/n.c*mean.c^2)))

#other effect size estimates
#library(compute.es)
#data <- data %>%
  #mutate(d=mes(mean.t, mean.c, sd.t, sd.c, n.t, n.c, level = 95, #cer = 0.2, dig = 2, , id = ID, data = data))
#check metafor

#data from ag and grazing studies that examined restoration in drylands
data <- read_csv("data/data.csv")
mydata<- data %>% 
 filter(disturbance == c("agriculture","grazing")) %>% 
 filter(!notes %in% "couldnt extract data")


mydata <- mydata %>%
  mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/n.t*mean.t^2) + (sd.c^2/n.c*mean.c^2))) #use only lrr now

Step 4. Summarize

Explore summary level data of all data. Explore aggregation levels that support the most reasonable data structure and minimize non-independence issues.

#evidence map####
require(maps)
world<-map_data("world")
map<-ggplot() + geom_polygon(data=world, fill="gray50", aes(x=long, y=lat, group=group))
#map + geom_point(data=paperdata, aes(x=long, y=lat)) +
  #labs(x = "longitude", y = "latitude") #render a literal map, i.e. evidence map, of where we study the niche in deserts globally

#add in levels and color code points on map####
#map of 178 articles
map + geom_point(data=data, aes(x=long, y=lat, color = paradigm)) + 
  scale_color_brewer(palette = "Paired") +
  labs(x = "longitude", y = "latitude", color = "")

#map of selected articles, agriculture and grazing disturbances
map + geom_point(data=mydata, aes(x=long, y=lat, color = paradigm)) + 
  scale_color_brewer(palette = "Paired") +
  labs(x = "longitude", y = "latitude", color = "")

#aggregation####
se <- function(x){
  sd(x)/sqrt(length(x))
}

data.simple <- mydata %>%
  group_by(study.ID, paradigm, technique, measure.success) %>%
  summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))

main.data <- mydata %>%
  group_by(study.ID, paradigm, intervention, outcome) %>%
  summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))


#EDA data####
simple.data <- mydata %>% group_by(study.ID, paradigm, technique, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
simple.data <- na.omit(simple.data)

parad.data <- mydata %>% group_by(study.ID, paradigm) %>% summarise(mean.rii = mean(rii), error = se(rii))
parad.data <- na.omit(parad.data)

tech.data <- mydata %>% group_by(study.ID, technique) %>% summarise(mean.rii = mean(rii), error = se(rii))
tech.data <- na.omit(tech.data)

success.data <- mydata %>% group_by(study.ID, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
success.data <- na.omit(success.data)


#active
active <- mydata %>%
  filter(paradigm == "active")

#viz for aggregation####
disturbance.data <- mydata %>% 
  group_by(study.ID,disturbance, paradigm) %>% 
  count()

disturbance.data2 <- disturbance.data %>% 
  group_by(disturbance,paradigm) %>% 
  count()

ggplot(na.omit(disturbance.data2), aes(disturbance,nn, fill=paradigm))+ 
  geom_bar(stat = "identity") + 
  coord_flip(ylim=0:44) + 
  scale_fill_brewer(palette = "Blues")

intervention.data <- active %>% 
  group_by(study.ID,intervention, paradigm) %>% 
  count()

intervention.data2 <- intervention.data %>% 
  group_by(intervention,paradigm) %>% 
  count()

#ggplot(na.omit(intervention.data2), aes(intervention,nn, fill=paradigm)) +
  #geom_bar(stat = "identity") + 
  #coord_flip(ylim=0:44) + 
  #scale_fill_brewer(palette = "Blues")

technique.data <- mydata %>% 
  group_by(study.ID,technique, paradigm) %>% 
  count()

technique.data2 <- technique.data %>% 
  group_by(technique,paradigm) %>% 
  count()

technique.data3 <- technique.data2 %>% 
  group_by(technique) %>% 
  count()


ggplot(na.omit(data.simple), aes(technique, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired")

ggplot(na.omit(data.simple), aes(measure.success, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired")

#better
ggplot(main.data, aes(intervention, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired") +
  labs(fill = "")

ggplot(main.data, aes(outcome, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired") +
  labs(fill = "")

Step 5a. EDA

Exploratory data analyses to understand data and QA/QC using Rii.

Step 5b. Models

Meta and conventional statistical models to explore relative efficacy.

Step 5. Synthesis stats

#p-value meta
#library(metap)
#all data for metas but cleaned for na's
#all data for meta cleaned
mdata.all <- mydata %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr))

#active only data for meta cleaned up
mdata <- mydata %>%
  filter(paradigm == "active") %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr))

#aggregated data for metas var estimated with central tendency
#note - could also bootstrap mean variance here instead of arithematic mean
simple.mdata <- mdata %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))

simple.mdata.2 <- mdata %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))  

simple.mdata.var <- mdata %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))

simple.mdata2.var <- mdata %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))

#metas with p-values####
#schweder(mdata$p)
#sumz(p, data = mdata)
#mdata %>%
  #split(.$intervention) %>%
  #purrr::map(~sumz(.$p)) 
#sumlog(mdata$p)

#metas with meta package on effect size measures####
library(meta)
#all data non-aggregated
#m <- metagen(lrr, var.es, studlab = ID, byvar = intervention, data = mdata)
#summary(m)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)

#Propose we only use random effect models because of the diversity of studies
m <- metagen(lrr, var.es, studlab = ID, byvar = paradigm, data = mdata.all)
summary(m) #active is positive and passive is negative so should not mix

## Number of studies combined: k = 496
## 
##                                          95%-CI             z p-value
## Fixed effect model   -0.1348 [-0.1348; -0.1348] -111701369.65       0
## Random effects model -0.0139 [-0.0855;  0.0578]         -0.38  0.7045
## 
## Quantifying heterogeneity:
## tau^2 = 0.2648; H = 12688836.57 [12688835.55; 12688837.59]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 12701660.72 [12701659.70; 12701661.74]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                     Q d.f. p-value
##  79698253852777280.00  495       0
## 
## Results for subgroups (fixed effect model):
##                      k                     95%-CI                    Q
## paradigm = active  333 -0.1348 [-0.1348; -0.1348] 79607575146295024.00
## paradigm = passive 163 -0.1574 [-0.1574; -0.1574]    90524299681441.44
##                     tau^2    I^2
## paradigm = active  0.2649 100.0%
## paradigm = passive 0.3666 100.0%
## 
## Test for subgroup differences (fixed effect model):
##                                   Q d.f. p-value
## Between groups      154406800798.70    1       0
## Within groups  79698099445976464.00  494       0
## 
## Results for subgroups (random effects model):
##                      k                     95%-CI                    Q
## paradigm = active  333  0.1947 [ 0.1033;  0.2862] 79607575146295024.00
## paradigm = passive 163 -0.3537 [-0.4882; -0.2192]    90524299681441.44
##                     tau^2    I^2
## paradigm = active  0.2649 100.0%
## paradigm = passive 0.3666 100.0%
## 
## Test for subgroup differences (random effects model):
##                      Q d.f.  p-value
## Between groups   43.68    1 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#should do a t-test right now of paradigm is diff from 0

#interventions
m <- metagen(lrr, var.es, studlab = ID, byvar = intervention, subset = paradigm == "active", data = mdata.all)
summary(m)

## Number of studies combined: k = 333
## 
##                                          95%-CI             z  p-value
## Fixed effect model   -0.1348 [-0.1348; -0.1348] -111668536.95        0
## Random effects model  0.1947 [ 0.1033;  0.2862]          4.17 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.2649; H = 15484891.12 [15484889.86; 15484892.37]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 15531744.14 [15531742.88; 15531745.40]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                     Q d.f. p-value
##  79607575146295024.00  332       0
## 
## Results for subgroups (fixed effect model):
##                             k                     95%-CI
## intervention = vegetation 221 -0.1348 [-0.1348; -0.1348]
## intervention = soil        75  0.7272 [ 0.7269;  0.7274]
## intervention = water       37  0.0000 [-0.0014;  0.0014]
##                                              Q  tau^2    I^2
## intervention = vegetation 79607575088050608.00 0.2649 100.0%
## intervention = soil                 2099026.81 0.0509 100.0%
## intervention = water                      2.34      0   0.0%
## 
## Test for subgroup differences (fixed effect model):
##                                   Q d.f. p-value
## Between groups          56145406.47    2       0
## Within groups  79607575090149632.00  330       0
## 
## Results for subgroups (random effects model):
##                             k                   95%-CI
## intervention = vegetation 221 0.1403 [ 0.0378; 0.2428]
## intervention = soil        75 0.4617 [ 0.3633; 0.5601]
## intervention = water       37 0.0000 [-0.0014; 0.0014]
##                                              Q  tau^2    I^2
## intervention = vegetation 79607575088050608.00 0.2649 100.0%
## intervention = soil                 2099026.81 0.0509 100.0%
## intervention = water                      2.34      0   0.0%
## 
## Test for subgroup differences (random effects model):
##                      Q d.f.  p-value
## Between groups   91.80    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

m <- metagen(lrr, var.es, studlab = ID, byvar = intervention, subset = paradigm == "passive", data = mdata.all)
summary(m)

## Number of studies combined: k = 163
## 
##                                          95%-CI           z  p-value
## Fixed effect model   -0.1574 [-0.1574; -0.1574] -2736465.45        0
## Random effects model -0.3537 [-0.4882; -0.2192]       -5.15 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.3666; H = 747523.89 [747522.42; 747525.37]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 272268.84 [272267.47; 272270.22]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                  Q d.f. p-value
##  90524299681441.44  162       0
## 
## Results for subgroups (fixed effect model):
##                             k                     95%-CI                 Q
## intervention = vegetation  54  0.2697 [ 0.2697;  0.2697]   775623742262.70
## intervention = grazing     12  0.4327 [ 0.4327;  0.4328]           6686.32
## intervention = soil        97 -0.7666 [-0.7666; -0.7666] 11085228014467.45
##                             tau^2    I^2
## intervention = vegetation  0.0075 100.0%
## intervention = grazing    <0.0001  99.8%
## intervention = soil        0.1148 100.0%
## 
## Test for subgroup differences (fixed effect model):
##                                Q d.f. p-value
## Between groups 78663447918024.95    2       0
## Within groups  11860851763416.47  160       0
## 
## Results for subgroups (random effects model):
##                             k                     95%-CI                 Q
## intervention = vegetation  54  0.0434 [-0.0088;  0.0956]   775623742262.70
## intervention = grazing     12  0.3481 [ 0.3412;  0.3550]           6686.32
## intervention = soil        97 -0.4835 [-0.5740; -0.3929] 11085228014467.45
##                             tau^2    I^2
## intervention = vegetation  0.0075 100.0%
## intervention = grazing    <0.0001  99.8%
## intervention = soil        0.1148 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   446.83    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#use random effects mean and var estimate to plot

#do t-tests here too

#outcomes
m <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "active", data = mdata.all)
summary(m)

## Number of studies combined: k = 333
## 
##                                          95%-CI             z  p-value
## Fixed effect model   -0.1348 [-0.1348; -0.1348] -111668536.95        0
## Random effects model  0.1947 [ 0.1033;  0.2862]          4.17 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.2649; H = 15484891.12 [15484889.86; 15484892.37]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 15555326.66 [15555325.39; 15555327.92]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                     Q d.f. p-value
##  79607575146295024.00  332       0
## 
## Results for subgroups (fixed effect model):
##                     k                     95%-CI                    Q
## outcome = soil    143 -0.1348 [-0.1348; -0.1348] 79607530173242944.00
## outcome = plants  108 -0.1990 [-0.1990; -0.1989]          68146638.29
## outcome = animal   11  1.0512 [ 0.9693;  1.1331]                 0.05
## outcome = habitat  71  0.5577 [ 0.5577;  0.5577]        3395032696.46
##                    tau^2    I^2
## outcome = soil    0.2649 100.0%
## outcome = plants  0.3814 100.0%
## outcome = animal       0   0.0%
## outcome = habitat 0.0739 100.0%
## 
## Test for subgroup differences (fixed effect model):
##                                   Q d.f. p-value
## Between groups       41509872733.84    3       0
## Within groups  79607533636422272.00  329       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI                    Q
## outcome = soil    143  0.1906 [ 0.0791;  0.3020] 79607530173242944.00
## outcome = plants  108  0.5161 [ 0.2783;  0.7538]          68146638.29
## outcome = animal   11  1.0512 [ 0.9693;  1.1331]                 0.05
## outcome = habitat  71 -0.4266 [-0.5776; -0.2755]        3395032696.46
##                    tau^2    I^2
## outcome = soil    0.2649 100.0%
## outcome = plants  0.3814 100.0%
## outcome = animal       0   0.0%
## outcome = habitat 0.0739 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   344.41    3 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

m <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "passive", data = mdata.all)
summary(m)

## Number of studies combined: k = 163
## 
##                                          95%-CI           z  p-value
## Fixed effect model   -0.1574 [-0.1574; -0.1574] -2736465.45        0
## Random effects model -0.3537 [-0.4882; -0.2192]       -5.15 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.3666; H = 747523.89 [747522.42; 747525.37]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 272268.84 [272267.47; 272270.22]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                  Q d.f. p-value
##  90524299681441.44  162       0
## 
## Results for subgroups (fixed effect model):
##                     k                     95%-CI                 Q   tau^2
## outcome = habitat  45  0.2697 [ 0.2697;  0.2697]   775623742237.89  0.0075
## outcome = plants   21  0.4327 [ 0.4327;  0.4328]           6703.49 <0.0001
## outcome = soil     97 -0.7666 [-0.7666; -0.7666] 11085228014467.45  0.1148
##                      I^2
## outcome = habitat 100.0%
## outcome = plants   99.7%
## outcome = soil    100.0%
## 
## Test for subgroup differences (fixed effect model):
##                                Q d.f. p-value
## Between groups 78663447918032.59    2       0
## Within groups  11860851763408.83  160       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI                 Q   tau^2
## outcome = habitat  45  0.0091 [-0.0455;  0.0637]   775623742237.89  0.0075
## outcome = plants   21  0.3483 [ 0.3415;  0.3552]           6703.49 <0.0001
## outcome = soil     97 -0.4835 [-0.5740; -0.3929] 11085228014467.45  0.1148
##                      I^2
## outcome = habitat 100.0%
## outcome = plants   99.7%
## outcome = soil    100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   464.58    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#t-tests for outcomes diff from 0 but you can see using the 95% CI what is different

#check metafor for interactions?? in these big models or are we ok??

#brainstorm on how to explore?? techniques


#save but cut all this.
#m.study <- metagen(lrr, var.es, studlab = study.ID, byvar = intervention, data = mdata)
#summary(m.study)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)

#aggregated data
#m1 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata)
#summary(m1)
#funnel(m1)
#metabias(m1)
#forest(m, sortvar = intervention)

#different var estimate
#m2 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata.var)
#summary(m2)
#funnel(m2)
#metabias(m2)
#forest(m, sortvar = intervention)


#m3 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata.2)
#summary(m3)
#funnel(m)
#radial(m3)
#metabias(m2)
#forest(m, sortvar = intervention)

#m4 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata2.var)
#summary(m4)
#funnel(m)
#radial(m4)
#metabias(m2)
#forest(m, sortvar = intervention)

#forest(m1, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m2, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m3, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m4, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#t-tests for lrr diff from 0
#mdata %>%
  #split(.$intervention) %>%
  #purrr::map(~sumz(.$lrr)) 

#effect sizes plots
#need better ci estimates
#ggplot(simple.mdata, aes(intervention, lrr)) +
 # ylim(c(-2,2)) +
#  geom_point(position = position_dodge(width = 0.5)) + 
 # labs(x = "", y = "lrr") +
 # coord_flip() +
 # geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = position_dodge(width = 0.5)) +
 # geom_hline(yintercept = 0, colour="grey", linetype = "longdash")

#ggplot(simple.mdata.2, aes(intervention, lrr, color = outcome)) +
 # ylim(c(-2,2)) +
  #geom_point(position = position_dodge(width = 0.5)) + 
 # labs(x = "", y = "lrr", color = "") +
 # coord_flip() +
  #geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = #position_dodge(width = 0.5)) +
  #geom_hline(yintercept = 0, colour="grey", linetype = "longdash")

#ggplot(mdata, aes(lrr, color = intervention)) +
  #geom_freqpoly(binwidth = 0.5, size = 2) + 
  #xlim(c(-10, 10)) +
  #labs(x = "lrr", y = "frequency", color = "") +
  #geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
  #scale_color_brewer()

#ggplot(mdata, aes(lrr, fill = intervention)) +
  #geom_dotplot(binwidth = 1) + 
 # xlim(c(-10, 10)) +
 # labs(x = "lrr", y = "frequency", fill = "") +
 # geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
 # scale_fill_brewer()